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Abstract 

We perform rescaled range analysis upon the signals measured by Dual Particle Dynamical 
Analyzer in gas-liquid two-phase turbulent jets. A novel rescaled range analysis is proposed 
to investigate these unevenly sampled signals. The Hurst exponents of velocity and other 
passive scalars in the bulk of spray are obtained to be 0.59±0.02 and the fractal dimension 
is hence 1.41 it 0.02, which are in remarkable agreement with and much more precise than 
previous results. These scaling exponents are found to be independent of the configuration 
and dimensions of the nozzle and the fluid flows. Therefore, such type of systems form a 
universality class with invariant scaling properties. 
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1 Introduction 



It is well known that picture of turbulent cascade, in which large eddies break 
down into smaller ones, is a multiplicative process. This hierarchical cascade in 
turbulence can be described by fractal geometry [Q], which is characterized quanti- 
tatively by the fractal dimension of self-similar structure of turbulence Both 
the experimental and theoretical aspects have been studied in past decades. 

Lovejoy investigated the fractal nature of satellite- and radar-determined cloud and 
rain areas covering 6 orders of magnitude of area sizes Js]]. The area-perimeter re- 
lation, introduced by Mandelbrot yQ, was found to hold with the fractal dimension 
Dj = 1.35 ±0.05. A theoretical analysis was proposed by Hentschel and Procaccia 
i^. They developed a theory of turbulent diffusion and obtained the natural conse- 
quence that 1.37 < Di < 1.41, which is in excellent agreement with the empirical 
results of Lovejoy @|. 



Another important experiment concerns the studies on the fractal facet of the turbulent- 
nonturbulent interface in turbulent flows by Sreenivasan and his coworkers ivilsllsl]. 
Prasad and Sreenivasan applied the laser-induced fluorescence technique to obtain 
the images of two-dimensional cuts of turbulent jets [lid. Illll. Applying the box- 
counting method, they estimated the fractal dimension of turbulence interface in 
the K range and found that D2 = 2.36 ± 0.05 for all fully turbulent flows Tak- 
ing into account the influence of local fluctuations in the Kolmogorov scale on the 
surface area due to the multifractal nature of the rate of dissipation, slight revision 
was made lll2n . In addition, Huang et al. proposed a modified box-counting method 
and found that D2 = 2 A3 ± 0.04 in the K range for round jets [|l3n . 



All the above experimental results are based on the box-counting method, which 
may lead to several disadvantages and difficulties. Difficulties associated with box- 
methods are typically attributed to the lack of a significant scaling range, low 
signal-to-noise ratio that limits reliable determination of level sets, and possible 
inadequate records resulting in poor statistics confidence. Several methods pro- 
posed to analyze box-counting statistics that do not presume power-law behaviors 
have been employed to address the first issue, which aim at determining the scaling 
range, if any, in an unbiased fashion. In addition, the linear correlation coefficient 
in log-log plots is not high enough, which results in relatively large standard devi- 
ations. It is well know that, for fractals with und erlying cascade process, there are 
logarithmic periodic oscillations lfl4l [l5. 16, 17, 18, 19]. The logarithmic periodic 
oscillations can be used to explain why different samplings lead to different esti- 
mates of fractal dimension. Zhou and Sornette proposed that canonical averaging 
of variety of samplings can be used to eliminate logarithmic periodic oscillations 
and get more reliable fractal dimension lilSll . 



Alternatively, rather than using the classic box-counting method, Zhou et al. per- 
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formed a rescaled range analysis upon signals measured by Dual PDA to determine 
the fractal dimension of coaxial turbulent jet, in which one-dimensional cuts are 
handled and the unequally spaced time series were pre-processed using averag- 



ing to form equidistantly spaced series 120(1 . The fractal dimension was found to 



be Di = 1.42 ± 0.07. Since the rescaled range analysis is conducted based on 
averaging different subseries at each scale, possible log-periodic oscillations are 
eliminated. However, the moving averaging interpolation approach is not satisfy- 
ing. In this paper, we shall develop a variant of the rescaled range analysis to make 
it suit directly for unequally spaced time series such that the corresponding results 
are more precise. 

This paper is organized as follows. In Sec. [2l we propose a generalization of the 
classical rescaled range analysis. Section [3] reports the experiential results and cor- 
responding fractal dimensions for a variety of experiments under different condi- 
tions. By comparing with interface dimensions of scalar field, physical interpreta- 
tion is presented in Sec. IH Section [5] concludes. 



2 Rescaled range analysis 



Rescaled range analysis, also termed as R/S analysis or Hurst analysis, was orig- 
inally developed by Hurst bill . This analysis is based on a new statistical devel- 
opment and provides an approach for analysis and characterization of time se- 
ries which has no underlying periodicity, yet retains long term correlation 1122 , 
23|, |24j, |25|l. The classic R/S analysis is performed on a discrete time series data 
set {ui : i = 1,2, ■ ■ ■ ,n} with the time tj uniformly- spaced, which is however not 
suitable for other time series that are not equidistantly sampled, such as data mea- 
sured by Dual PDA. We thus generalize the classic R/S analysis from equidistant 
sampling to uneven sampling. 



Define a continuous function y : [a, b] ^ 3?. In a certain sense, graph 

graph (y) = {{t,y (t)) : t e [a,b]} 



(1) 



may be regarded as a 2-dimensional fractal in the plane {t,y). Falconer has pre- 
sented a method to estimate the fractal dimension of functional graph [)26.1. Let 
[ti, T2] C [a, b]. The time span s = r2 — n of ?/ is referred to as "lag" [1231 . 12411 . Then 
the mean of y on [ti, T2] is 



{y) 



1 



T2 



T2 — Ti Jti 

Define the cumulative deviation to be 

rt 



y (t) dt. 



X{t)= [y (t) - {y)] dt 



(2) 



(3) 



Tl 
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and the cumulative range to be 



R= sup X{t)- inf X{t). (4) 



The standard deviation is 



S={^- f[y{t)-{y)fdiY. (5) 

I T2 — n Jti ) 

Therefore, as a generalization of the R/ S analysis for discrete time series, we as- 
sume that R/S is self-affine and scales with respect to s as a power law 

R/S cx {t2 - nf = s", (6) 

where H is termed as the Hurst exponent. 

Partition the interval [ri, T2\ into n divisions with ti = < ti < ■ ■ ■ < tn = T2. 
Then Vi = y (ti) is a unequally spaced time series. Using rectangular approximation 
for integration, one can descretize Eqs. dH)-® as foUows.The mean of y can be 
regarded as a time-weighted average 

1 " 

{y) = _ & X {ti - , (7) 

and the cumulative range becomes 

i?= max{X,}- mm {X,}, (8) 

where 

i=l 

is the cumulative deviation. Taking into account the weighted standard deviation 

^ = {^-TT^E (y^ - {y)f (t^ - ^-1)} ' ' (10) 

we have 

R/ScK{tn-tof = s''. (11) 

If ti — = At is a constant, the cumulative deviation Xj given in ^ is At 
times that in the conventional R/S analysis for uniformly spaced time series data. 
However, it has no impact on the estimation of Hurst exponent using Eq. (fTT)) . The 
Hurst exponent H of the time series is then evaluated from the slope of the straight 
line fitted to these points. The estimated Hurst exponent is related to its fractal 
dimension D by the relation [|3|] 

D = 2-H. (12) 
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The R/ S analysis is a simple and robust method for fast fractal estimate with as 
few as 30 data points Mm- 



There are a variety of algorithms for performing R/ S analysis. Consider a subseries 
{Ui, Vi+i, ■ ■ ■ , Ui+s-i}, where i is the departure of the sub-series and s is the length 
of the sub-series with l^i<i + s — \ ^n. Different selection of i and s 
results in different algorithm. In this paper, we propose to adopt a random selection 
algorithm, in which m starting points i are chosen randomly based on a uniform 
distribution on the interval [1, n — s + 1]. An average of the m values of R/ S with 
the same lag s is calculated, which is denoted by R{s)/S (s). The Hurst exponent 
can be calculated according to Eq. (fTTl) . 



3 Experimental measurement of fractal dimensions 

3.1 Experimental set-up 

The Dual PDA is based on a novel concept that yields a higher measurement accu- 
racy and performs non-intrusive measurements of the velocity, diameter and transit 
time of spherical particles, droplets and bubbles suspended in gaseous or liquid 
flows, particularly for spray analysis and other investigations of liquid atomization. 
The underlying principle of phase Doppler anemometry is based on light-scattering 
interferometry and therefore requires no calibration. The measurement point is de- 
fined by the intersection of two or three pairs of focused laser beams and the mea- 
surements are performed on single particles as they move through the measurement 
volume. Particles thereby scatter light from the laser beams, generating an optical 
interference pattern. A receiving optics placed at a well-chosen off-axis location 
projects a portion of the scattered light onto multiple detectors. Each detector con- 
verts the optical signal into a Doppler burst with a frequency linearly proportional 
to the particle velocity. The phase shift between the Doppler signals from different 
detectors is a direct measure of the particle diameter. We can obtain the arrival time 
and transit time simultaneously. 

The Dual PDA receiving probe contains four receiving apertures integrated into 
one single optical unit. The Dual PDA detector configuration (2 standard and 2 
planar), combined with sophisticated validation routines, is not susceptible to un- 
wanted effects resulting from the Gaussian light intensity distribution in the mea- 
surement volume. Misinterpreted size measurements due to trajectory effects are 
therefore eliminated. The front optics module for 3D-PDA configurations simpli- 
fies the alignment procedure considerably. Screwed onto the PDA receiving probe 
and connected to the transmitter optics by a dual-fibre link, the front optics module 
generates the third pair of laser beams with adjustable beam intersection and focus. 
The received signals are fed to one of Dantec's advanced signal processors, which 
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delivers results to a PC. The instrument chart of the Dual PDA system is shown in 
Fig.[TJ The power of the laser generator is 2W. The focal length of the transmitting 
and receiving lenses are both 500mm. 
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Fig. 1 . Experimental set-up of the Dual PDA system 



3.2 Experimental conditions 



In the experiments, we used tri-passage coaxial nozzles. The liquid phase water 
moves through the mid-passage, while the gaseous phase, say air, N2 and CO2, 
passes through the inner and outer passages. In the experiments, the gaseous medium 
is changeable and the fluid flow rates of the two phases are adjustable. The mea- 
surement points distribute through out the spray zone. For a fixed experimental 
configuration and given fluid flow, we record the arrival time t, transit time T, axial 
velocity U, radial velocity V and drop diameter d of 20000 drops moving through 
the measurement point. In experiment (a), the flows of nitrogen in the inner and 
outer passages are respectively 0.5 and 3.0m^/h, while those in the experiment (b) 
are 0.69 and 3.9m^/h, respectively. 

For a fixed measurement point at (0, 0, —26), the data rate is 20.8KHz, the mean ve- 
locity is 33.81m/s, and the r.m.s. velocity is 9.51 m/s. The Kolmogorov microscale 
7] is calculated from the signals according to 



1/4 



IhUdU/dtf 



(13) 



where Ui is the mean speed at the measurement point, u is the kinematic viscosity 
of air, and U is the streamwise component of velocity. The resulting value of rj is 
0.013cm. The Taylor microscale A = 0.41cm is calculated according to 



A= . (14) 

[{dU/dtt 
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where U' is the r.m.s. of velocity fluctuations. Hence, the Taylor-scale Reynolds 
number is 

Re;, = Ui\/u = 273 , (15) 

which is moderate. 

The power spectrum of the velocity signal is shown in Fig. [2] which is obtained by 
the Dual PDA Processor using the Gabor Fourier transform. A power-law scaling 
E{k) ~ k"f^ with an exponent /? close to 5/3 is observed over a substantial range of 
about 0.7 decades. Since the sampling frequencies are not high enough, the highest 
frequency is in the inertial range, which implies that we are dealing with the K- 
range in this work. 




1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 

Fig. 2. Energy spectrum of velocity signals 
3.3 Hurst exponents and fractal dimensions 

First of all, we analyzed the velocity signals arising from experiments by random 
selection algorithm. A typical diagram of ln(i?/S') versus ln(s) is shown in Fig.[3l 
There is a transient regime for ln(s) < —6, which exhibits a clear dropdown com- 
pared with the scaling regime. The slope of the trend line in the scaling regime 
gives the Hurst exponent H = 0.5926. When adopting the approximation that the 
measurement is evenly spaced, the Hurst exponent estimated in the same scaling 
range is H = 0.5518, which is significantly smaller than the real value of H. This 
phenomenon is systematically observed for other experimental realizations. 

According to the trend of Hurst exponents from the experiments, we can roughly 
classify the turbulent spray zone into two parts: the transient region and the steady 
region. In addition, the span of the transient zone is found to be proportional to the 
spray speed. Two typical diagrams of variation of Hurst exponent along the axial 
direction are shown in Fig. SI It is obvious that there exists transient region near 
the outlet of the nozzle. The span of transient zone of (a) is about 2 ~ 3mm, while 
that of (b) is about 3 ~ 4mm. One can see that H decreases along the streamwise 
direction near the outlet of the nozzle, which shows that the signals close to the 
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Scaling regime 



-10 -8 -6 -4-2 2 

In (5) 

Fig. 3. A typical diagram of Iti{R/S) versus ln(s) 

nozzle have stronger persistency. This is induced by the stronger interactions among 
drops moving through the measurement points. After a short transient period, the 
four types of H arrive a "steady" state, which fluctuates somewhat randomly near 
the average 0.59. We should point out that the measurement range in the experiment 
with a distance from the outlet of the nozzle of about 200mm is much wider than 
what we have presented in Fig. lU 
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Fig. 4. Typical diagrams of Hurst exponent evolution along axial direction 

The Hurst exponents of axial and radial velocities are plotted in Fig. \5\ The exper- 
iment conditions and/or the spatial position of measurement are different to each 
other. The hurst exponent is independent of experiment condition and the spatial 
position. It is easy to find that, there are several points with relatively higher H of 
the axial velocity. These signals are nonstationary and are not available for R/S 
analysis, which will be addressed later. Thus we withdraw these points. It follows 
that 

i/u = 0.60 ± 0.02 (16) 

and 

ifv = 0.59 ± 0.02. (17) 

In fact, the existence of nonstationary signals does not affect the mean values, but 
they will increase the standard deviations. 
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Fig. 5. Hurst exponents of velocity signals 

Similarly, the Hurst exponents of transit times and drop diameters are shown in 
Fig. [6l Nonstationary signals, especially of transit time, appear again. One can also 
find that, when investigating a signal with high H in Figl6l those signals from the 
same measurement record of the investigated signal have high H as well. With- 
drawal of these nonstationary signals follows that 



Ht = 0.60 ± 0.02 



(18) 



and 



Ha = 0.59 ± 0.02 



(19) 
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Fig. 6. Hurst exponents of transit time and drop diameter signals 

It is obvious that the Hurst exponents of different signals are identical to each other. 
The difference among the four types of H is within the scope of experiment error. 
Consequently, the considered signals are also self-affine among each other. The 
histogram of the Hurst exponent distribution of all signals is shown in Fig. U\ The 
majority of Hurst exponents concentrate around the mean value 0.59. The averaged 
Hurst exponent can thus be calculated as 



if = 0.59 ± 0.02 . 



(20) 
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We can hence obtain the fractal dimension of the signals that 

D = 1.41 ±0.02 . 



(21) 




Fig. 7. Histogram of Hurst exponents of all signals 

In a nutshell, the Hurst exponent and the corresponding fractal dimension of the 
signals are independent of the spatial position of the measurement position. That 
is, the investigated turbulent jet is fractally homogenous in the main bulk of the 
spray zone. Moreover, we can say that such two-phase flows form a universal class 
with a universal Hurst exponent, since the fractal dimensions are invariant with 



the changes of nozzle configuration, fluid medium, and flow rates as well [|20n . 
However, we should point out that the Hurst exponents change for other type of 
two-phase flows. For example, coherent structure appears when measuring signals 
of hydrogen jet into air [i28l l29i1 . Such a self-organized structure is expected to 
strengthen the long-term dependence and decrease the fractal dimensions. 



4 Discussions 



4. 1 Nonstationary signals 



One may take it for granted that the computed Hurst exponent of a fixed record 
is not reproducible when a repeated calculation is carried out. However, we would 
like to point out the stability and repeatability of the random sub-series algorithm, 
which has been verified by repeating computations such that the resulting Hurst 
exponents change very slightly compared with the previously calculated values. 
Thus, it seems that the surprising high values of Hurst exponents in Fig.|5]and Fig.|6] 
are unavoidable. Fortunately, we found that signals corresponding to "unexpected" 
high values of H are nonstationary which can be used to distinguish these signals 
from the rest. This is why we have got rid of experimental points with high H. 
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A random process, say the axial velocity U(t), is called stationary if identical rules 
generate the process U (t) itself and all the processes deduced from U (t) by a time 
shift, namely, all the processes of the form U{t + s).ln experiments, we have to 
control the conditions to be fixed through out the measurement of one signal. Con- 
sequently, the mean of axial velocity component, denoted as {U (s)) where the 
averaging is performed over the time and s is the sub-series length, must be a con- 
stant. For nonstationary signals in the experiments, an obvious change of {U (s)) 
is detected in the plot of {U (s)) against arrival time. In the present case, the non- 
stationary signal can be regarded as a patchwork of several stationary ones with 
different means. We find that, such a scrambling leads to an upward in the tail of 
large lag s, namely, the resulting Hurst exponent becomes higher. From a mathe- 
matical point of view, (f[^^ for these nonstationary signals is dependent of t and 



thus R/S analysis is not applicable pSi] 



It is the usual situation that the flow flux decreases or increase suddenly. Therefore, 
in order to verify the validity of the claim, we performed numerical computations 
upon artificial series u (ti) which are generated using 

, U (ti) , if 1 ^ z ^ 10000 

u{U) = { (22) 

U (ti) + An, if 10001 ^ I ^ 20000 

from the real series U (U), where An is an adjustable non-zero parameter. Denote 
if" for the Hurst exponents of the artificial series. We find that H"^ (An) is approxi- 
matively equal to (—An) in a statistical sense, and H"^ increases with increasing 
I An I in the sampling time intervals for all series. We use the same signal analyzed 
in Fig.[3]to show the effect of non stead by comparing the results of raw signal with 
manipulated signals (An ^ 0). For the sake of better presentation, we treat all the 
signals as evenly spaced. The results are shown in Fig. [H The circles represent the 
real series, while the squares and pluses correspond to An = 10 and An = —10, 
respectively. The solid line is used to guide the eyes. We can see that the real series 
has scaling behavior over the whole sampling interval. Nevertheless, the artificial 
series have a more narrow scaling range and an upward end, which results in the 
increase of the Hurst exponent. 

Imagine a sudden drop in the flux of nitrogen in the outer passage such that {U (s)) 
decreases. (T(s)) will increase considerably due to the strong dependence of tran- 
sit time on the axial velocity, while {d{s)) and {V{s)) increases slightly. These 
changes result in the in- step increases of H. 



4.2 Comparison with interface dimension 



It is natural to think that the passive scalars, the transit time and drop size, are 
strongly affected by the turbulent velocity field. Since the passive scalars are de- 
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Fig. 8. Comparison of rescaled range analysis on the real time series (o) and manipulated 
signals with Au = 10 (□) and Au = -10 (+) 

pendent of the turbulent field, the fluctuation of passive scalars shows the situation 
of turbulent field. If the investigated carriers are dense enough such that they form 
a continuous turbulent-nonturbulent interface, we can conjecture that the fractal 
dimension of the turbulent-nonturbulent interface depend strongly upon the fluctu- 
ation of the velocity and the passive scalars. We will find that the fractal dimensions 
of interfaces and signals are essentially identical. 

It is obvious that, if the drops in the spray zone is dense enough, the spray can be 
regarded as a "cloud". Furthermore, if the disperse phase in the spray aggregates to 
form larger bulk, interfaces appear. Such crossovers illustrate the intrinsic relation- 
ship among these systems. It is convenient to assume that there are local imaginary 
interfaces in the investigated system. It is well known that the area S of the fractal 
interface can be expressed as yO 

S^r]'^~^\ (23) 



We write C for a scalar difference characteristic of large scales, such as the r.m.s 
value, and 7] for the Kolmogorov scale. It is shown that the characteristic scalar 
gradient across interfaces is of the order C/?7 dst]. According to the K41 theory, we 
have 

ri ~ (24) 
Therefore, the flux of momentum across the interface can be written as 

flux = uSC/t] ~ Z/(7-3i32)/4_ (25) 

It is now well known that the growth rate of turbulent flows of a given configu- 
ration are independent of fluid viscosity, which is referred to as Reynolds number 
similarity ^ . It follows that 

D2 = 7/3. (26) 

There are two possible corrections to D2 according to Eq.(l25l). By taking into ac- 
count fluctuations in r] due to the multifractality of the rate of energy dissipation. 
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one may do such correction by computing the mean value of C/rj. Further correc- 
tion is made by considering the estimate of interface area 5" due to fluctuations in 
rj. The corresponding results is in essential agreement with the experimental out- 
comes. 

Comparison of our results with the previously mentioned situations follows that the 
fractal dimension D of the signals is related to the interface dimension D2 by 

D = D2-1. (27) 

In other words, D can characterize the turbulent-nonturbulent interface as well as 
D2. Actually, the signals fluctuate weaker at the beginning of the spray and the 
interface, if exists, is much smoother. This can account for the increasing trend of D 
in the transient region near the nozzle. In addition, the intrinsic relevancy between 
the fluctuations of signals and the interface can be also explained by Taylor's frozen 
flow hypothesis. Certainly, the fluctuations of D close to the mean is due to the 
variation of mean velocity within a single record. 



5 Conclusions 

In this paper, we perform rescaled range analysis on the signals measured by the 
Dual PDA in gas-liquid two-phase turbulent jets. We generalize the classical R/S 
analysis to continuous form and then discretize it to make it suitable for unequally 
spaced time series. The fractal dimensions of the signals of axial and radial ve- 
locities, transit time and drop diameter under different experimental conditions are 
obtained both in the transient region and the bulk of the spray zone. The fractal 
dimensions of signals of different physical quantities are identical to each other. 

In the transient region, since the liquid phase is being sped up by the high-speed 
gaseous phase, the fluctuations become more and more remarkable along the spray 
direction. As a consequence, the fractal dimensions increase streamwisely. The 
length of the transient region decreases with increasing gas flow rate. In the bulk of 
the spray region, we have D = 1.41±0.02 for all signals and experiment conditions, 
which is in excellent agreement with interface fractal dimension by experiments 
and theoretical derivations and is not effected by the drop size distribution. Since 
the fractal dimension is invariant for experimental configurations and conditions, 
the gas-liquid two-phase turbulent jets investigated in this work form a universality 
class with an invariant exponent. 

Certainly, there are still open problems. To address why the Hurst exponents (and 
corresponding fractal dimensions) of different variables at different measurement 
points are identical, we have to uncover the underlying physics. However, we can 
only provide a qualitative explanation rather than quantitative expressions. But to 
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make further progress, we would need some kind of theory or model to give the 
power law relationship between the different variables. In addition, it is necessary 
and interesting to test the log-periodic oscillations in the turbulent signals, which 
may provide an important step towards a direct demonstration of the Kolmogorov 
cascade or at least its hierarchical imprint [Il8l.ll9n. 
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